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Solving the Schrodinger eigenvalue problem by the imaginary time 
propagation technique using splitting methods with complex coefficients 
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The Schrodinger eigenvalue problem is solved with the imaginary time propagation technique. The separability 
of the Hamiltonian makes the problem suitable for the application of splitting methods. High order fractional 
time steps of order greater than two necessarily have negative steps and can not be used for this class of 
diffusive problems. However, there exist methods which use fractional complex time steps with positive real 
parts which can be used with only a moderate increase in the computational cost. We analyze the performance 
of this class of schemes and propose new methods which outperform the existing ones in most cases. On the 
other hand, if the gradient of the potential is available, methods up to fourth-order with real and positive 
coefficients exist. We also explore this case and propose new methods as well as sixth-order methods with 
complex coefficients. In particular, highly optimized sixth order schemes for near integrable systems using 
positive real part complex coefficients with and without modified potentials are presented. A time-stepping 
variable order algorithm is proposed and numerical results show the enhanced efficiency of the new methods. 
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I. INTRODUCTION 



We consider the eigenvalue problem for the stationary 
Schrodinger equation (SE) subject to asymptotic bound- 
ary conditions (h = m = I), 



where 



H<f>i(x) = Ei<j)i(x), 



H = T + V{x) 



0,1,2, 



-A + V(x) 



(1) 



(2) 



V(x) denotes the interaction potential and A is the 
Laplacian operator. Since the Hamiltonian if is a Her- 
mitian operator, its eigenvalues E± are real valued, and 
its corresponding real eigenfunctions <f>i(x) form a basis 
of the underlying Hilbert space. This particular prob- 
lem has attracted great interest among theorists and 
practitioners^— due to its relevance for the understand- 
ing of the atomic and molecular structure of matter. 

A widely used approach to solve this problem is based 
on using the corresponding time-dependent Schrodinger 
equation in imaginary time (t = — ir), whose formal solu- 
tion is given by the evolution operator cxp(— tH). In this 
way, in general, any initial condition, under the action of 
cxp(— tH), converges asymptotically to the ground state 
solution when t — > oo. Notice that the evolution opera- 
tor exp(— tH ) has the same eigenfunctions as the prob- 
lem fl}-©. This technique is usually referred to as the 
imaginary time propagation method (ITP for short). In 
this setting, only the action of exp(— tH) on a wave func- 
tion has to be computed^. 



Since the operators e~ rV and e _rT can be exactly com- 
puted in the coordinate and momentum space, respec- 
tively, the operator splitting technique involving a com- 
position of these exponential operators with appropriate 
coefficients can be used to approximate e~ rH . The com- 
putational cost depends on the number of changes be- 
tween these coordinates which are cheaply performed by 
Fast Fourier transforms (FFT). 

The ITP method can be regarded as an analog of the 
well-known power method in numerical linear algebra^. 
In this sense, one may also consider the inverse power 
method: instead of the iterative application of the ex- 
ponential operator exp(— tH), the scheme v n +\ = (H — 
Ei)~ 1 v n , n = 0, 1, 2, . . . is used. This iteration is known 
to converge after normalization to the eigenvector with 
eigenvalue closest to Ei. Although faster convergence 
than for the ITP method can be observed for an accu- 
rate initial guess Ei m Ei, in general, the algorithm needs 
more iterations until convergence^. 

One could also consider polynomial approximations to 
cT tB , such as those provided by the Chebyshev and 
Lanczos methods^. In that case, one has to evaluate 
products of H with vectors, which can also be computed 
using FFTs. However, the number of products (or the 
degree of the polynomials used) for the whole time inte- 
gration, i.e., the computational cost of the method, has 
to be proportional to \\H\\, and this norm depends on the 
number of grid points. Splitting methods on the other 
hand do not have such a restriction, their computational 
cost is simply proportional to the number of time-steps 
n, and this is irrespective of the number of mesh points 
since c~ hv and c~ hT arc well defined operators for all 
h > 0. Even though the norm of T increases for a finer 



mesh the value of h does not have to be reduced to make 
the algorithm stable. As a result, splitting methods scale 
better than polynomial approximations like Chebyshev 
or Lanczos, in agreement with the results observed in 
the numerical experiments ir£. 

The operator splitting technique has also some limita- 
tions. In particular, splitting methods of order p > 2 re- 
quire negative time-step a 10 ' 11 and the instabilities caused 
thereof are analogous to the ones for the integration of 
a diffusion equation backwards in time. If it is feasi- 
ble to compute the gradient of the potential V, general- 
ized splitting methods allow to build methods with posi- 
tive coefficients up to fourth orde r 12 ' 13 , but higher order 
methods also use negative time-steps. In this paper we 
propose methods to overcome the order barriers for both 
cases by using complex time-steps. Splitting methods 
can be tailored to particular equations to achieve bet- 
ter performances and we present criteria based on near- 
integrability that apply to a wide range of problems and 
thus yield highly efficient high order schemes. The ob- 
tained methods outperform the existing splitting schemes 
when high accuracy is desired and could be appropriate 
for elaborating a variable order algorithm. We also report 
some numerical experiments illustrating the efficiency of 
the new methods. 



II. THE IMAGINARY TIME INTEGRATION METHOD 
FOR THE SCHRODINGER EQUATION 

An important property of the Hcrmitian operator H 
is that (choosing properly the origin of the potential) its 
eigenvalues < Eq < E\ < . . . are real and nonnegative, 
and the corresponding cigenfunctions fa can be chosen 
to form a real orthonormal basis on its domain. The 
problem $Q originates from the time-dependent SE 

ig- t ^(x, t) = H^(x, t), i>{x, 0) = i/j (x). (3) 



A Wick rotation of the time coordinate, t 
forms ((3]) into a diffusion type equation 



-ir, trans- 



d_ 



if>(x,r) = Hi/j(jc,t), ^(»,0)=Vo(aO, (4) 



with formal solution ip(x,r) 



-tH 



%jj{x,Q). After ex- 



panding the initial condition ipo in the basis of eigen- 
functions fa, 

ipo(x) = J^c 4 fa(x), Ci = (fa(x) ip(x,0)) , 

i 

where (• | •) is the usual L 2 scalar product, the time evo- 
lution of (Ql is given by 

i>(x, r) = e~ rH ^(x, 0)=J2 e ~ TEx c * <M X )- ( 5 ) 



Asymptotically, for a sufficiently long time integration, 
we get %\)(x,t) — > e~ rE ° cofa since the other exponen- 
tials decay more rapidly. The convergence rate depends 



of course on the separation of the eigenvalues. For sim- 
plicity, we restrict ourselves to the non-degenerate case 
Eq < Ei. If there is degeneracy, it converges to a linear 
combination of eigenfunctions, and repeating this process 
with different initial conditions one can obtain a complete 
set of independent vectors of the subspace which can be 
orthonormalized. 

Normalization of the asymptotic value yields the eigen- 
function fa and the corresponding eigenvalue is com- 
puted via Eq = (fa\Hfa). Excited states can be ob- 
tained by propagating different wave functions simulta- 
neously (or successively) in time and using, for example, 
the Gram-Schmidt orthonormalization or diagonalizing 
the overlap matrix-^. 

For simplicity in the presentation, the spatial dimen- 
sion is set to one unless it is explicitly stated, but our 
results also apply to higher dimensions. 

The problem is further simplified by assuming x G [a, b] 
with the interval [a, b] sufficiently large such that the 
wave function and all its derivatives of interest vanish 
at the boundaries. For numerical computations, the in- 
finite dimensional domain of H has to be truncated, 
which is done by discretizing the spatial coordinate x: 
wc fix N equally spaced grid points x% = xq + kAx, k = 
0, 1, 2, . . . , N — 1, with a — xq and b = x^. In this 
way, the interval is divided into N subintervals of size 
Ax=(b- a)/N. 

The potential V is represented in this grid by a diag- 
onal matrix and the periodicity of the system ("0(a) = 
-0(6) = 0) allows for the use of spectral methods (in 
space) for the calculation of T, namely the Fast Fourier 
Transform after which the matrix representation of T also 
becomes diagonal. The computational costs for the appli- 
cation of V and T to a vector are thus proportional to N 
and NlogN operations, respectively. In a d-dimensional 
space with N mesh points on each dimension, their costs 
are proportional to N d and N d log N, respectively. 



III. SPLITTING METHODS FOR THE SCHRODINGER 
EQUATION 

To approximate the time evolution ([S]), i.e., the com- 
putation of cT tH acting on a vector, we propose to use 
compositions of the operators e~ rV and e~ rT evaluated 
at different times. A first example is provided by the 
well-known Strang splitting 



*M = e-? v e- hT e-? v , 



(6) 



verifying ^ 



[2] _ n -hH 



0(h 3 ) with h = At. Higher 



order approximations can be obtained by a more general 
composition 



* 



W 



n 



e -a,ihT e -bihV 



(7) 



where "J 



[p] 



-liH 



0(h p+1 ) if the coefficients a,, bi are 
chosen such that they satisfy a number of order condi- 



tions (with m sufficiently large). It is well-known, how- 
ever, that methods of order greater than two (p > 2) nec- 
essarily have negative coefficient ;} 10 ' 11 ' 14 (a simple proof 
can be found ini£). While this is usually not a prob- 
lem for the coefficients bi, having negative Oj coeffi- 
cients makes the algorithm badly conditioned (in the 
limit TV — > oo). Recall that this restriction applies to 
a more general class of equations: whenever the operator 
only creates a semigroup of solutions, the integration has 
to be computed forward in time. 

Composition methods with coefficients bi positive are 
also convenient for unbounded potentials, e.g., V{x) = 
x , since negative values of bi can generate large round- 
off errors in the exponential e~ biV at the boundaries if 
the interval-size of the spatial discretization is not appro- 
priately chosen and the potential takes exceedingly large 
values. 

Splitting methods are particularly appropriate for the 
numerical integration of this problem since the choice of 
the time step, h, is not affected by the mesh size. Taking 
a finer mesh (i.e., a larger value of N) does not necessar- 
ily lead to a smaller time step, and the extra computa- 
tional effort originates only from the FFTs, whose cost is 
N \og(N) (or N d log( N ) in a d-dimensional problem with 
N points on each coordinate). In contrast, the Cheby- 
shev polynomial method approximates the exponential 
as 



-tH 



-P(H/\\H\\ 



Y,ak(p)Pk(rH/p), (8) 



fc=0 



where p = r||if||, afc(p) = 2Ik(p), fc > and ao(p) = 
Io (p) , with /„ being the modified Bessel functions of the 
first kind and Pk is the Chebyshev polynomial of degree 
k. The truncation index m has to be such that m > p to 
reach accuracy. If we take into account that 



\H\\ 



V,„ 



V m 






n 
Ax, 



N 2 



we observe that the number of FFTs is proportional 
to TV 2 , and then its computational cost scales as TV 2 x 
N d log(N). Similar considerations apply to other polyno- 
mial approximations such as Taylor or Lanczos methods. 
In fact, in the numerical experiments with d = 3 reported 
in£, splitting methods scale as TV 33 while Lanczos meth- 
ods scale as N 5 - 3 for a range of values of N, in agreement 
with the previous estimates. 

Nevertheless, from the results in£ one might conclude 
that these polynomial approximations scale better than 
splitting methods for computing a larger part of the spec- 
trum. Then, splitting schemes are the methods of choice 
when only relatively few eigenvalues are desired, espe- 
cially if they are needed with high accuracy (implying a 
finer mesh, i.e., a larger number of grid points). With this 
aim in mind, we next explore several families of splitting 
methods to get approximations to ([5]). 

One possible approach to derive the order conditions to 
be satisfied by the coefficients aj, bi consists in applying 



the Baker-Campbcll-Hausdorff formula to the composi- 
tion (|7||, which we assume consistent (J2i a i = J2i b i = 
1)±£. Thus we get *Jf ] = exp(-hU), with 



H=T + V + hf 2il [T,V] 
+ h 2 (f 3A [T,[T,V}]- 



f 3 ,2[V,[T,V]]) 



(9) 



where /).., arc polynomials of degree i in the coefficients 
dk , bk and the symbol [T, V] stands for the commutator 
of the operators T and V. Condition /2.1 = leads to 
second order methods, and this can always be achieved by 
taking a left-right symmetric composition in ([7} because 
all even terms automatically vanish. Methods of higher 
orders require in addition /31 = /3 a = 0. Taking into 
account consistency, these equations can be written asii 



/3,i : £ cubjdk = -, 



l<i<j<fe<m 



/s,2 : £ hdjbk = -. 



(10) 



(11) 



l<i<j<k<m+l 



These two conditions imply that at least one of the a^ as 
well as one of the bi become negative (seei£ and references 
therein), so that only methods of order two can be used 
for this problem. 

There are several possibilities to circumvent this limi- 
tation, and in the following, we enumerate some of them. 

a. Modified potentials. If the kinetic energy opera- 
tor in (J4|) is quadratic in momenta, then the nested com- 
mutator 



[V,[T,n = (W(x)f(W{x)) 



(12) 



is diagonal in coordinate space. For this reason, (TT2"]) 
is usually called modifying potential. In consequence, 
[V, [V, [T,V]}} = and we can replace the terms e - b > hv 
in by the more general operator e ~ b i h v-cih [V,[t,v\] 
involving two parameters. As a result, condition (|11[) 
becomes 



l<i<j<k<rn+l »=1 



(13) 



This equation can always be satisfied with a proper 
choice of the coefficients Cj, so that the constraints on the 
coefficients a,,6i reduce to the single condition /3.1 = 0, 
allowing for positive coefficients. In addition, solutions 
with positive q coefficients also exist. A first example is 
the 4th-order compositio n 12 ' 18 



*M= e -^V 



™ V -%[V,[T,V]] , 



(14) 



It turns out, however, that 6th-order methods using the 
operator (|12p necessarily have some negative coefficients 



.19 



b. N ear- integr able systems. When the Hamiltonian 
can be considered as a perturbed system, i.e., H = H + 
eV E {x) with an exactly solvable part H = T + Vo{x) 
and a small perturbation eV £ (x), it is advantageous to 
split the Hamiltonian into the dominant part Hq and its 
perturbation eV e . For example, if one is interested in the 
lower excited states, which evolve near the minimum of 
the potential, it can be useful to separate the quadratic 
part and to treat the remainder as a perturbation since 
the harmonic oscillator has a simple and fast solution 
using FFTs^. 

Notice that in this case, the commutator 



[eV e ,[H ,eV e ]} 



(W e (x)) T (W E (*)) 



depends only on the coordinates and modified potentials 
can also be applied as before. Then, all compositions 
remain the same except for replacing T by Hq and V by 
eV £ . 

With the additional information that one part of the 
operator is significantly smaller than the other, it is clear 
that the error expansion for a consistent method ^h can 
be asymptotically expressed as 



*/, - c 



-hH 



= ^2 Yl et > k £lhk > as C 1 ' e ) ~* (°> °)> 



analytic scmigroup a 23 ' 24 . It is worth noticing that the 
propagator exp(zA) (z e C) associated with the Lapla- 
cian is well defined (in a reasonable distributional sense) 
if and only if Re(z) > 022. More generally, it is possible 
to extend the semigroup related with parabolic partial 
differential equations into a sector in the right half plane 
ofC2£. 

Many of the existing splitting methods with com- 
plex coefficients have been constructed by applying the 
composition technique to the symmetric second-order 
leapfrog scheme (JB|). For example, a fourth-order inte- 
grator can be obtained with the symmetric composition 



*W = ^ ^[2] ^[2] 



where 



1 



2 - 2 1 /3 e 2«^/3 ' 



/? 



2l/3g2ifc7r/3 
2 _ 2l/3 e 2'*-^/3 ' 



(15) 



(16) 



and k = 1,2. In both cases, one has Re(a),Re(/3) > 0. 
Higher order composition methods with complex coeffi- 
cients and positive real part can be found ir*22r— , where 
several numerical examples are also reported. 



i>l k>s t 



where the Si start from the first non-vanishing error co- 
efficient e Su k- We say that ^h is of generalized order 
(si, S2, ■ ■ ■ , s m ) (where s\ > S2 > ■ ■ ■ > s m ) if the local 
error satisfies that 



*; 



-hH 



0(eh Sl 



81 + 1 



2j,«a+l 



e z h 



■e m h Sm+1 ). 



Thus, for a method *f? h 
error reads 



of generalized order (8, 2), the 



* 



-hH 



iA sh 9 + e 3 , 2 e 2 h 3 + 0(e 3 h 3 



This class of schemes can also be applied in several other 
situations. For instance, suppose one takes a sufficiently 
fine mesh. Then ||T|| ^> ||V|| and the previous consid- 
erations apply (with Ho = T). Also, if V(x) = x n for 
a sufficiently large n, then the virial theorem (cf> \T\ </>) = 
((f>\VV(x)x\(t>) lcads'to (T) =n(V). 

c. Complex coefficients. A third possibility consists 
of considering complex coefficients in the composition ([?P 
(with or without modified potentials). In other prob- 
lems where the presence of negative real coefficients is 
unacceptable, e.g., in those concerned with differential 
equations defined in semigroups, cf. Q, the use of high 
order splitting methods with complex coefficients hav- 
ing positive real part has shown to possess some ad- 
vantages. In recent years a systematic search for new 
methods with complex coefficients has been carried out 
and the resulting schemes have been tested in different 
settings: Hamiltonian systems in celestial mechanics^, 
the time-dependent Schrodinger equation in quantum 
mechanics^ and also in the more abstract setting of evo- 
lution equations with unbounded operators generating 



IV. NEW SPLITTING METHODS FOR THE ITP 
PROBLEM 

In this section, we carry out a systematic search of 
methods within the classes (a)-(c) above enumerated. We 
only consider symmetric methods and, since T and V 
have qualitatively different properties, we analyze both 
TVT- and VTV-type compositions, defined as 

V^b] _ c -a t hT e -bihV c -a 2 hT . _ . e ~a 2 hT Q -bihV g-oihT 

(17) 



and 



M _ -bihV - ai hT -b 2 hV . , . e -b 2 hV e -a 1 hT e -b 1 hV 



*)f = e-" inv e~ ainl e 



(18) 



respectively. In principle, both compositions have the 
same computational cost for the same number of expo- 
nentials. Nevertheless, due to a projection step to the 
real part after each full time-step, only in the VTV com- 
position we can concatenate the last map in the current 
step with the first stage in the next one. The TVT com- 
positions thus require two additional FFTs in comparison 
with the VTV composition, and this is accounted for in 
the numerical experiments. 

The methods we obtain are classified into two fami- 
lies: (I) methods without modified potentials and (II) 
methods with modified potentials. For each class we dis- 
tinguish between methods for general problems (with the 
unique constraint that [V, [V, [T, V]]] = 0) and methods 
for near-intcgrable problems (when the main dominant 
part contains the kinetic energy). 



We have explored both TVT and VTV compositions 
with different number of stages. In some cases we con- 
sider extra stages to have free parameters for optimiza- 
tion. When the number and complexity of the order con- 
ditions is relatively low, we get all solutions. We then 
select the solutions having all of their coefficients with 
positive real part. Finally, we choose the solution which 
minimizes 



X!(i°' 



(19) 



and/or minimizes the absolute value of the real part 
of the coefficients appearing at the leading error terms. 
These methods are subsequently tested on several nu- 
merical examples. After this process, we collect a num- 
ber of schemes offering the best performance for most of 
the problems considered. In practice, however, one has 
to bear in mind that the relative performance between 
different methods depends eventually on the particular 
problem considered, the desired accuracy, the initial con- 
ditions, etc. 



A. Methods without modified potentials 

TVT and VTV compositions with 3 up to 9 stages 
have been analyzed. To simplify the notation, we denote 
compositions (|T7|) and ([15]) as 



Tn m = 0,1 bi a 2 
Vn m = b\ ax b 2 



a 2 b 1 a 1 , 
b 2 oi &i 



respectively. Here n indicates the order (or generalized 
order) of the method and m corresponds to the num- 
ber of stages, i.e., the number of bi coefficients in the 
TVT composition or the number of a» coefficients in the 
VTV composition. The coefficients of the selected TVT 
methods are collected in Table U whereas those corre- 
sponding to the TVT methods are displayed in Table [Til 
In all cases, the three coefficients in the middle of each 
composition (in bold) are not included since they can be 
obtained from the previous ones by consistency (all of 
them have also positive real part). The coefficients are 
expressed with 18 digits for simplicity. If required, the 
same tables collecting the coefficients with 25 digits are 
available from the authors upon request. 



where the Eij are chosen to form a basis of the algebra 
L l (T, V) of commutators of length i. The chosen basis 
elements relevant for our exposition are 



£■3,1 = 


= [T, [T,V]], 


Ez. 2 - 


= [V,[T,V]}, 


■E-5,1 = 


= [T,[T,[T,[T,V]]]], 


E§. 2 - 


= [V,[T,[T,[T,V\]]], 


^5,3 = 


= -[T,[V,[T,[T,V]]]], 


E5.4 - 


= [V[V,[T,[T,V]]]], 


■EV,1 = 


--[T,[T,E 5!l ]}, 


E7.2 - 


= \V,[T, Es.i]]- 



Here we summarize some of the methods which have been 
analyzed: 

a. 3-stage compositions. A 3-stage composition has 
sufficient parameters to build 4th-order methods. There 
is one real solution and two complex solutions (conjugate 
to each other). For example, the VTV method corre- 
sponds to the composition (fT5J) when W h * is given by ((6]). 
The TVT version is obtained by interchanging T and V . 

b. 5-stage compositions. Fourth-order methods 
with two free parameters can be obtained using 5-stage 
symmetric compositions. These two parameters can 
be used to build methods of effective order 6 (i.e., 
4th-order methods that arc conjugate to 6th-order 
methods by a near-identity change of variables). This 
requires to impose some additional constraints on the 
leading error terms, /sj, j = 1,2,3,4. Specifically, 
these are f^\ — /g^ = and /s^ + ^4 = 02£. We have 
found six solutions for the TVT composition and three 
solutions for the VTV composition with coefficients 
having positive real part. The solutions with smaller 
error terms at order 5 are denoted by T45 and VI5. 

c. 7-stage compositions. In principle, there are suf- 
ficient parameters to build 6th-order methods with 7 
stages. For the TVT composition there are 11 solutions 
with all coefficients having positive real parts. The solu- 
tion leading to a minimum value of the norm of the error 
at order 7 is denoted by T67 in Table |TJ 

With respect to the VTV composition, the best 
method we have found is given in Table HTl V67. It is iden- 
tical to the most efficient sixth-order method obtained by 
Chambers^, where it has been presented as a symmetric 
composition similar to (fT5|) but with 7 stages instead of 

3, and with ^ h given by ([6]). 



2. Methods for near-integrable problems 



1. Methods for general problems 

Analogously to (|9]), the symmetric compositions (fl7|) 
and p8p can be formally expressed as a single exponen- 
tial &£' = cxp(-hTL) with polynomials fcj in a*;, bi mul- 
tiplying commutators Eij: 

00 dim(L'(T,V)) 
n = T + V+ J2 Jl h'^fijE^, 

i— p+1 j — 1 



Proceeding analogously as before, we arrive at the fol- 
lowing methods. We recall that in all compositions one 
should replace T by Hq and V by eV e . 

a. n-stage compositions of generalized order (2n,2). 
This class of compositions has real and positive 
coefficients^. A 4-stagc VTV composition of general- 
ized order (8, 2) is given by scheme V84M^ R in Table HVl 
with ci = 0. 

b. 5-stage compositions. To build a method of gen- 
eralized order (8,4) the following conditions must be sat- 
isfied by a consistent and symmetric method: /31 = 



i*3,2 = /s,i = /7.1 = 0. It requires at least 5 stages, 
and in this case only one solution with all coefficients 
having positive real part is found both for the TVT and 
VTV compositions. The coefficients of these methods, 
denoted by T84s and V84s, are collected in Table U and 
Table HH respectively. 

c. (8,6,4) methods. To build a (8,6,4) method, the 
coefficients of a consistent and symmetric method must 
satisfy the following order conditions: /^i = fop = 
/s,i = ha = /s,3 = fr,i = 0. They therefore require 
at least 7 stages. In this case, it is possible to get all so- 
lutions. Scheme T8647 corresponds to the solution mini- 
mizing (|19|) . whereas V864y provides the minimum value 

of 1/5,3 + hA- 

d. (8,6) methods. Increasing the number of stages 
to 9 we have two free parameters, which are used to sat- 
isfy in addition the following conditions: /s : 4 = /^ = 0. 
In this way, methods of generalized order (8,6) and effec- 
tive order (10,8,6) are obtained. Two efficient schemes 
correspond to T869 and V869 in Table [J and Table HO 
respectively 2 ^. 



B. Methods with modified potentials 

Fourth-order methods incorporating modified poten- 
tials do exist with real and positive coefficients. In fact, 
2- and 3-stage schemes have been extensively studied 
(se o 13 i 19 ). Methods of generalized order (n, 4) also ex- 
ist with positive real coefficients 2 ^. Here we construct 
new methods of generalized order (6,4) and (8,4) with 
this property and generalize the treatment to 6th-order 
schemes with complex coefficients. In all cases, we take 
compositions TVT and VTV with up to 5 stages and 
denote them as 

TnM m = a\ {b 1 ci)a 2 ■ ■ ■ a 2 (61 ci)a\, 
VnM m = (61 ci) 01 {b 2 c 2 ) • • • (62 c 2 ) a\ {h a). 

Here, the parenthesis is used to help counting of the num- 
ber of exponentials, and the letter M indicates that the 
methods use modified potentials. Notice that the num- 
ber of free parameters can differ for the TVT and VTV 
sequences with the same number of exponentials because 
the exponent of a modified potential contains two pa- 
rameters. The coefficients of the selected methods are 
collected in Table [TTT] and Table [TV] for the TVT and 
VTV compositions, respectively. Again, the three coeffi- 
cients in the middle of each composition which appear in 
bold (only for the coefficients a^ and bi) are not included 
since they can be obtained from the previous ones by 
consistency (all of them have positive real part). 



1. Methods for general problems 

a. 4-stage compositions. Under the restriction of 
having real positive coefficients, we have obtained the 



fourth-order VTV method OMF-4M, already discovered 



,13 



(eq. (36) therein). 



The VTV composition allows one to build 6th-ordcr 
methods, whereas the TVT needs an extra stage. There 
is only one solution (and its complex conjugate) with all 
coefficients having positive real part. It is denoted by 
V6M4 in Table HVl 



2. Methods for near-integrable problems. 

We first consider (n, 4) methods with real and positive 
coefficients. For schemes of generalized order (8,6) we 
collect only complex solutions with positive real part. 

a. (6,4) methods They require at least 3 stages to 
satisfy the following order conditions: /3 : i = $■$$ = 
/51 = 0. The coefficients a\ and bi correspond to the 
methods (6,2) obtained ir>2i (without modified poten- 
tials). We have also considered methods with 4 stages 
in order to have additional free parameters. As pre- 
viously mentioned, there is the same number of order 
conditions as parameters to get a method of order 6 for 
the VTV sequence, but there are no solutions with co- 
efficients being real and positive. To get a sixth-order 
method the following conditions must also to be satis- 
fied: /5,2 = /5,3 = /s,4 = 0- The coefficients Cj only 
appear in /s^ and /s.4 and can only be used to cancel 
these terms. The VTV sequence has three free parame- 
ters which can be used to annihilate /s^ and ^4 and to 
minimize the absolute value of /s^ under the constraint 
that all coefficients must be real and positive. The TVT 
sequence has only two free parameters which can be used 
to annihilate f$^ and to minimize the absolute value of 
the dominant term, f$.2, under the same constraint on 
the coefficients. The best methods we have obtained are 
denoted by T64M 4 and V64M 4 . 

b. (8,4) methods They require at least 4 stages. The 
coefficients a% and bi correspond to the methods (8,2) 
without using modified potentials and obtained in21. 
There is one coefficients c, in the TVT composition which 
can be used to cancel /5,3, and two coefficients c, in the 
VTV composition which can be used to annihilate /^ 
and /5.4. The solution with c 2 = 03 = was already 
obtained in2£. We have collected the corresponding coef- 
ficients for this method, V84M^ R , in Table EI We have 
also considered methods with 5 stages in order to have an 
additional free parameters. There is the same number of 
order conditions as parameters to get a method of order 
(8,6) (which would be of order 6 for a general problem) 
but, obviously, there are no solutions with coefficients 
real and positive. As in the previous case, the term f^, 2 
can not be zeroed using real positive coefficients. Then 
in both TVT and VTV compositions we have chosen the 
method which, while having real and positive coefficients, 
minimize its absolute value. The best methods we have 
obtained arc denoted by T84M 5 and V84M 5 . 

c. (8,6) methods They require at least 5 stages and 
do not admit real and positive solutions for the coeffi- 



TABLE I. Compositions TVT without modified potentials. 

T4 5 = Oi bi Q2 b 2 a 3 b 3 a 3 b 2 01 &i Oi 

al = 0.087177040194768004 + 0.069259433759657453; 
61 = 0.174354080389536008 + 0.138518867519314907; 
a 2 = 0.208397663533315769 + 0.052478646956409802; 
b 2 = 0.242441246677095530 - 0.033561573606495303; 

T67 = ai fei a 2 62 03 b 3 a4 \>4 a.4 63 a3 62 a 2 bi ai 



TABLE II. Compositions VTV without modified potentials. 



ai = 0.060011447243322255 

61 = 0.115650213164896285 
a 2 = 0.113925528178022672 

62 = 0.128728705301158783 
a 3 = 0.166160237564322898 

63 = 0.186892287438669413 



+ 0.020040690498800685i 
+ 0.044802659670983506; 

- 0.037781051835206749i 

- 0.123466712987851772i 

- 0.065264806973654210; 
+ 0.001153105240840699; 



T84 5 = ai 61 a 2 b 2 a 3 b 3 a 3 b 2 a 2 61 ai 

al = 0.071401131540044698 + 0.010155431019886789; 
61 = 0.178696854264631978 + 0.028197506313218021; 
a 2 = 0.236383805190074736 + 0.070427007139534522; 
b 2 = 0.198453474708154649 + 0.082962314733854963; 

T8647 = ai bi a 2 b 2 03 63 &4, h^ a.4 63 03 b 2 a 2 bi a\ 



ax = 0.055705821110864236 + 
61 = 0.115779449626990422 + 
a 2 = 0.118843282163492564 - 
b 2 = 0.129128920804026450 - 
a 3 = 0.158591515575195578 - 
b 3 = 0.184643464154438944 - 



0.018670384565085049; 
0.046131356173382847; 
0.024151805322796634; 
0.119039413303774209J 
0.076302551893579599; 
0.003053761445376182; 



T869 = ai fei a 2 b 2 a 3 b 3 04 &4 a 5 b 5 a 5 64 04 b 3 a 3 b 2 a 2 61 a\ 



ai = 0.042257897299860339 
61 = 0.094894869367770736 
a 2 = 0.095260398471830494 
b 2 = 0.097374660381711248 
a 3 = 0.099960578944766657 

63 = 0.118584793520055816 
a 4 = 0.148695530402608487 

64 = 0.136865119760326031 



- 0.014215780224181831; 

- 0.037963806472588094; 
+ 0.004518725891475591; 
+ 0.088518877931710497; 
+ 0.090271995071312563; 
+ 0.038356250608401259; 
+ 0.011438117187614089; 

- 0.023587404969570006; 



V4 5 



61 ai b 2 a 2 b 3 a 3 b 3 a 2 b 2 ai 61 



61 = 0.087177040194768004 + 0.069259433759657453; 
ax = 0.174354080389536008+0.138518867519314907; 
b 2 = 0.208397663533315769 + 0.052478646956409802; 
g 2 = 0.242441246677095530 - 0.033561573606495303; 

V67 = 61 ai b 2 a 2 &3 a 3 b4 a4 b4 a 3 b 3 a 2 b 2 ai bi 



61 = 0.058450018777330642 
01 = 0.116900037554661284 
b 2 = 0.123229569418374773 
a 2 = 0.129559101282088262 
63 = 0.158045797047111040 
a 3 = 0.186532492812133817 



+ 0.021714127308030170; 
+ 0.043428254616060341; 

- 0.040280678786016125; 

- 0.123989612188092593; 

- 0.060441090739009958; 
+ 0.003107430710072675; 



V84 5 = 61 a\ b 2 a 2 b 3 a 3 b 3 a 2 b 2 a\ 61 

61 = 0.052472525516129026 - 0.010958940842458138; 
01 = 0.175962140656732362 - 0.054483056228160557; 

62 = 0.246023563332753880 - 0.125228547924834352; 
a 2 = 0.181259898687454283 - 0.034864508232090522; 

V8647 = b\ ai b 2 a 2 b 3 a 3 \>4 8L4 b4 a 3 b 3 a 2 b 2 ai 61 



61 = 0.060017770752528926 - 
ai = 0.108904710931114447 - 
b 2 = 0.067017987316853817 + 
a 2 = 0.106594114300156182 + 
b 3 = 0.189300872388005476 + 
a 3 = 0.204897016414416105 + 



0.009696150746907738; 
0.075700232434276860; 
0.003927567742822542; 
0.139651903644940761; 
0.091055103879530385; 
0.009719057955143112; 



V869 = 61 ai 62 a 2 b 3 a 3 64 04 b 5 a 5 b 5 04 64 03 63 a 2 b 2 ai b\ 



61 = 0.032497706037458608 
ai = 0.087895680441261752 
b 2 = 0.094180923422602148 
a 2 = 0.095351855399045611 

63 = 0.101132953097231180 
a 3 = 0.121865575594908413 

64 = 0.160941382119434892 
a 4 = 0.141506882718462097 



+ 0.010641310380458924; 
+ 0.036052576182866484; 
+ 0.023866875362648754; 

- 0.065128376035135147; 

- 0.112201757337044841; 

- 0.054974002471495827; 

- 0.016127643896952891; 
+ 0.024607229046524026; 



cients and we are forced to consider complex solutions. 
We have found only one solution with positive real part 
in the coefficients for both TVT and VTV compositions. 
The coefficients for the methods denoted by T86M5 and 
V86M5 are given in Table IIIII and Table IIV1 respectively. 



V. NUMERICAL EXAMPLES 

A. Efficiency of the methods 

As test bench for the numerical methods, we con- 
sider in the following two qualitatively different cases, the 
Poschl- Teller potential and a perturbed harmonic oscilla- 
tor, the latter being a classic example of a near-integrable 
system and of practical interests^. These two problems 
can be numerically integrated using modified potentials. 
However, we compare the relative performance of the 



methods (with and without modified potentials) sepa- 
rately in order to study the performance of the methods 
when it is not feasible to compute the gradient of the 
potential. 

The numerical integration proceeds as follows: starting 
from random initial data, we iterate with fixed time-step 
until the sufficiently large final time T = 100 and com- 
pare the result with the exact solution, ip(T), which has 
been obtained by integrating with a much smaller time 
step. The spatial interval is fixed for all experiments to 
[—10,10] and is discretized with N = 512 equidistant 
mesh points. At each step, we project the obtained vec- 
tor to its real part and normalize it to one in £2^), i-e., 
given the method w£ and initial conditions, u n € M. N , 
we compute u„+i as 

u n+1 = ^ u n ; 
then, since w„+i is a complex vector (but 0(h p ) away 



TABLE III. Compositions TVT with modified potentials. 



T64M4 = oi (bi ci) a 2 (b 2 c 2 ) a 3 (b 2 c 2 ) a 2 (bi ci) «i 



ai = 0.088344743225421001 
61 = 0.214113867876328652, 
a 2 = 0.274436837849719332, 



a = 0.000457570336450475 
c 2 = 0.000967171438516093 



T84M 5 = 01 (&i ci) a 2 (62 c 2 ) a 3 (b 3 c 3 ) a 3 (b 2 c 2 ) a 2 (61 ci) ai 



01 = 0.058520963359694865 
61 = 0.145381537601615725, 
a 2 = 0.207903047442871771 
b 2 = 0.244351408696638327, 
c 3 = 0.000938105701711153 



Ci 



0.000245906549261228 



c 2 = 0.000259178561419125 



T86M5 = 01 (61 ci) a 2 (62 c 2 ) a 3 (b 3 c 3 ) a 3 (b 2 c 2 ) a 2 (h ci) ai 



01 = 0.063556051997493102 
bi = 0.156939525347224563 
ci = 0.000133739181746125 
a 2 = 0.208998817231756322 
b 2 = 0.222383136675982213 
c 2 = 0.000484323504408882 
c 3 = 0.000179180363327321 



+ 0.010606890396680920i 
+ 0.027931306200415819i 
+ 0.000085540153220213i 
+ 0.040240203826523395i 
+ 0.026033262090035938i 
+ 0.000241671051573332i 
- 0.00085830441303451 \i 



TABLE IV. Compositions VTV with modified potentials. 

V6M4 = (bi ci) Oi (b 2 c 2 ) a 2 (b 3 c 3 ) a 2 (b 2 eg) Oi (bi ci) 

bi = 0.070024595602928186 + 0.015533573983094094i 
a =0 

01 = 0.210073786808784559 + 0.046600721949282283i 
b 2 = 0.280190696984660355 + 0.045035902822440560i 
c 2 = 0.000906238223230531 + 0.000614631850819718i 
c 3 = 0.000526292468950591 - 0.001566397557699757J 

V64M 4 = (61 ci) ai (b 2 eg) a 2 (b 3 c 3 ) a 2 (b 2 c 2 ) Oi (bi Ci) 



bi = 0.064342100381112671, 
01 = 0.204839710862229619 
b 2 = 0.293632462299624640, 
c 3 = 0.001525544947900195 



ci = 0.000217958066128410 
c 2 = 0.000601872205939880 



from a real vector) we project on the real space by re- 
moving the imaginary part 

u n +i = Re(u n+ i) 

and finally we normalize the solution 



«n+l 



where the norm is given by 



V84M 5 



bi = 0.042308451243127365, 
ai = 0.142939324267716184 
b 2 = 0.219303568753387110, 
a 2 = 0.242474508234531493, 



Ci 



0.000232966269565498 



c 2 = 5.56677120231130 ■ 10" 
c 3 = 0.000794490777479431 



V84M4- R = (bi ci) ai (b 2 c 2 ) a 2 (b 3 c 3 ) a 2 (b 2 c 2 ) ai (bi ci) 



61 = 1/20, 
b 2 = 49/180, 



f-i 



_ 3 861-79lV21 
129600 ' 

c 2 = 0, c 3 = 



Ol 



1/2- ^3/28 



bi = 0.046213625838152095 - 0. 
ci = 0.000035830461339520 + 0. 
01 = 0.152650950104799817 - 
b 2 = 0.224258052678856384 - 0. 
c 2 = 0.000338053435041382 - 0. 
a 2 = 0.226364275186039762 - 
c 3 = 0.000408311644874003 + 0. 



007824529355983108i 
00007437085768542 U 
,030279967163699065i 
050879282402761772i 
000490508913279372J 
,016537249619936515i 
000484371967433683i 



u n+ i 



11 



n+ll 



\w\f 



Ax 



N-l 

3=0 



(w ,...,w N -i) e 



»A' 



We take as the computational cost the number of Fourier 
transforms necessary until the final time. In addition, 
the methods using complex coefficients are penalized by a 
factor 2 in the computational cost, which comes from the 
use of complex Fourier transforms instead of real FFT. 
We repeat the numerical integrations for different values 
of the time step, i.e., h = T/M for different values of M. 
We take as the approximate solution, <fi(T) = u n in each 
case and measure the error in position as 

error = \\<tp(T) - <j>(T)\\. 

This procedure will allow us to determine the efficiency 
of the new splitting methods, which will depend on the 
desired accuracy, and thereby choose the methods which 
are most appropriate for implementation with a more ef- 
ficient algorithm that is based on variable time step and 
order. We distinguish two types of problems, on the 
one hand, methods that include modified potentials, the 
reference methods being Chin-4M (JT4J), OMF-4MH and 
V84M^ R ^ given in Table ED and on the other hand, 
methods without modifying potentials with the reference 
methods V822£, the fourth order complex triple-jump 
scheme (|15p . referenced as Yoshida 4 and a 6th-order 
complex coefficient method by Chambers^. 



(61 ci) ai (b 2 c 2 ) a 2 (b 3 c 3 ) a 3 (b 3 c 3 ) a 2 (b 2 c 2 ) a\ (h ci) 



Poschl-Teller potential 

We have chosen the well-known one-dimensional 
Poschl-Teller potential for the availability of analytic so- 
lutions of the eigenstates and also for its spectrum that 
does not only consist of eigenvalues, 



V86M5 = (bi ci) 01 (b 2 eg) «2 (b 3 e 3 ) a 3 (b 3 eg) «2 (b 2 eg) 01 (bi ci) 



H 



1 d 2 
'2'dx 1 



A(A + 1) 



'sech(x) 



1 



(20) 



with A(A + 1) = 10. The results of our computation are 
shown in Figure Hal The higher order of the complex 
coefficient methods outweighs their extra cost starting 
from moderate accuracy. The optimizations of the error 
terms can be clearly appreciated in the comparison with 
the 4th order triple-jump (|15j). When we consider meth- 
ods with modified potentials, we observe that the new 
methods show only slight improvements with respect to 



the method OMF-4M since both parts of the splitting T 
and V are of comparable size. As the desired precision 
is increased, the new sixth order methods dominate in 
efficiency. 



2. Perturbed harmonic oscillator 

To illustrate the benefits of methods designed for near 
integrable systems, we use the Hamiltonian 

19 2 1 92 



1 a 2 



and split in a large part Huo = ^^dx 1 ~^~ 2 an d a 

small part eV^cc). The trap frequency is set to u = 1 
and the perturbation el^ is given by the Poschl- Teller 
potential in (|20]). with A(A + 1) = 2/5. The harmonic 
part -Hho can be solved exactly via an exact splitting 
using Fourier transforms, cfi^, where it is shown that 

— i5H HO = g -i Tf tan(^)a; 2 — i^j sin(<5w) p 2 — i§ tan(^r)x 2 

~2 

for |<5o;| < 7r and p = — -gs- After the substitution 
6 = —ih, we have 

— /iZ/ho -— ^- tanh(^)a — ^j sinh(/iu) p — ^tanh(4p)a; 

for |Im(/i)w < 7T and Rc(/i) > (for numerical stabil- 
ity) and the perturbation part is easily propagated after 
discretization by the exponential of a diagonal matrix. 
In this setting, the higher order in the small parameter 
is amplified and the efficiency plots in Figure llbl indi- 
cate that the new methods outperform the existing ones 
when high precision is sought and overall when modi- 
fied potentials are allowed. We observe in both examples 
that, when modified potentials can be computed with- 
out exceedingly large computational cost, they should be 
used. 

Further numerical experiments show that the efficiency 
curves are independent of the mesh size, i.e., the norm of 
T, and the cost only increases as Alog(iV) as expected. 
The reason for this can be understood by following the 
evolution of the state vector along the iterations of the al- 
gorithm. Whereas in the beginning one has a non-smooth 
configuration Uq, after a few steps the vector Ui is close 
to an eigenstate and thus smoothened. 

It is important to remark that the methods proposed in 
this work can be implemented in an algorithm which uses 
variable step, variable order, variable mesh size and vari- 
able simple-double precision. The best implementation 
can depend on the class of problems to be solved. For 
illustration, we present an implementation with variable 
time steps. 



B. Variable step method 

The previous examples show that for low accuracies 
and large time steps, the (8,2) method (with real co- 



efficients) performs best. However, if we allow for vari- 
able time steps, as proposed in£»£, the computational cost 
is drastically reduced. We propose an improved time- 
stepping algorithm that is based on two different estima- 
tors for the eigenvalue. 

Recall that fixing the time-step and iterating to con- 
vergence will yield an eigenvector with the error being of 
the order of the method 0(h p ) since we are computing 
exactly the spectrum of a perturbed Hamiltonian. As- 
sume now that we are close to convergence, i.e, one has 
obtained an eigenvector u n = vq + O (h p ) and we consider 
the decomposition in the basis of exact eigenvectors Vi of 
H, 



N-l 



i n = 2J (kvi, where 



N-l 

E 

i=0 



\di? = l. 



It is clear that di — 0(h p ), i > 1 and due to the nor- 
malization do = l + 0(h 2p ^). Then, an energy estimation 
is given by 

E K1 = u T n Hu n = E + 0{h 2p ). 

Alternatively, the energy can be estimated by the loss of 
norm in each time step, 

u n+ i = e- hH u n + <D{h p+1 ) = e- hE °v + 0(h p+1 ), 

and then 

Eh a = log (H*"+ 1 ID = Eo + chP + O (ft** 1 ) . 
h 

Combining both expressions yields an error estimate for 
the energy, 

AE h = E K2 - E h , x = ch p + 0(h p+1 ). 

The convergence in energy is measured by comparison 
with the previous time step, 

SEJl = E% A - E^- 1 = dh 2p + 0(h 2p+1 ). 

The time stepper then works as follows: starting from 
a large step size, the time step is decreased by a factor 
1/2 whenever the actual reduction in energy of the itera- 
tion 8E falls below the the maximally reachable precision 
AE, i.e., \SE\ < (AE) 2 and the iteration is terminated 
once the error estimate AE has reached a given tolerance. 
This algorithm is summarized in Table [V] 

For the numerical experiments, we use the same config- 
urations as for constant time step, however, the iteration 
is terminated when convergence in energy is reached at 
AE < 10~ 10 . The iterations are initialized with ran- 
dom normalized data and a time step of r = 10. The 
results are displayed in Figure [5a] for the Poschl- Teller 
potential and in Figure I2bl for the perturbed harmonic 
oscillator with the same parameters as in the fixed-step 
size experiments. The error is measured as the £2 norm of 
the difference between the current value of the algorithm 
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(a) Unperturbed Poschl- Teller potential 



(b) HO perturbed by a Poschl- Teller potential 



FIG. 1. (color online) In the first row, efficiency curves (error vs. number of FFTs) for methods without force evaluations 
are presented, with the new methods (triangles) performing best for high accuracies. The middle rows depicts methods based 
on modified potentials. In the left column, T86M5 intersects with V84M 5 at precision 10 -13 , whereas it already takes place 
at 10 -9 for the right column. In the bottom row, the random initial conditions (green), the ground states (black) and the 
potentials (dashed blue), scaled by 1/5 to fit the axis, are shown. 



ip(t) and the exact ground state 4>(T) as in the previous 
experiments, error = \\ip(t) — <fi(T)\\. 

As expected, it is apparent that lower order methods 
show better smoothing behavior for the first steps, when 
the wave function is still rough (recall that the algorithm 
is initialized by a worst-case wave function). For higher 
precisions, the new methods clearly outperform the ex- 
isting ones, with the sole exception of the unperturbed 
setting with modified potentials, where the globally op- 
timized OMF-4M method can hardly be improved un- 
less extremely high precision is sought and the 6th or- 



der methods of Table IIVI and IIIII become favorable (not 
shown) . 

The results indicate that for low precision, i.e., for the 
first iterations, a lower order method should be used and 
then, after a certain precision is reached, e.g., when the 
higher order methods exhibit their superiority the algo- 
rithm should change to the optimal method, either V8647 
or V86M5 until convergence. Further preliminary exper- 
iments on this adaptive order strategy have shown that 
there is plenty of room for optimization, e.g., by changing 
the initial step-size, adjusting the step-size by a differ- 
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(a) (color online) Unperturbed case 



(b) (color online) Perturbed case 



FIG. 2. Evolution of precision in the li norm of the position vector with the variable time step algorithm described in Sec. 
IV Bl As in Fig. [TJ the top row gives the results for standard methods whereas the bottom rows shows methods with modifying 
potentials. 



Time step controller 

Input: h, wo, tolerance, n := 0; 
do 

n := n + 1; 

compute : u n := Re(vt '£ u n -i); 

normalize : u n :— u n /\\u n \\; 

compute : 5E 

if (\SE\ < (AE) 2 ) then h ~ h/2 
while \AE\ < tolerance 
Output: eigenvector u n 



and AE 



For excited states, one expects an even better perfor- 
mance of the new methods since several states have to 
be computed to high precision in order to avoid error ac- 
cumulation and the gains of the new methods are thus 
amplified. We have confirmed this conjecture by numer- 
ical experiments. The results thereof are omitted in the 
manuscript since they do not contribute insight beyond 
the presented experiments: they are qualitatively identi- 
cal. 



TABLE V. Adaptive time stepping algorithm for a given v| _ CONCLUSIONS AND OUTLOOK 



W 



method ^ 



ent factor or by modifying the control criterion. Each of 
which has certain advantages and disadvantages, depend- 
ing on the initial conditions and the range of precision. 



We have studied the Schrodinger eigenvalue problem 
by the imaginary time propagation method and proposed 
splitting schemes with complex coefficients that can over- 
come the order barrier for parabolic problems since the 
coefficients have only positive real parts. The obtained 
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sixth order methods are clearly superior to any classical 
ones for high precisions. On the other hand, when the 
gradient of the potential can be cheaply evaluated, the 
high order methods with complex coefficients are efficient 
only at very high accuracies due to the double cost caused 
by complex arithmetic. 

We have proposed different high order methods to 
reach highly accurate results. An efficient implementa- 
tion should take into account, for example, a preliminary 
time integration on a coarse mesh using simple precision 
arithmetic in order to get, as fast as possible, a smooth 
and relatively accurate solution from a random initial 
guess, and next consider a refined mesh using arithmetic 
in double precision. For simple precision arithmetic and 
low accuracies, it suffices to consider only low order meth- 
ods, and when higher accuracies are desired we turn to 
double precision, variable time step and variable order 
methods. The best algorithm could depend on the class 
of problems to solve. 

It is also important to remark that the form of the ex- 
ponent allows that the techniques presented in this work 
can also be transferred to other areas whenever splitting 
is appropriate and the integration has to be performed 
forward in time, e.g., statistical mechanics of quantum 
systems, where one has to compute the Boltzman oper- 
ator exp(— j3H), with j3 = (kT)^ 1 or quantum Monte- 
Carlo simulations^. 
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